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p  We  derive  new  upwind  type  finite  difference  approximations  to  systems  of 
nonlinear  hyperbolic  conservation  laws.  The  general  technique  is  exemplified 
by  the  potential  flow  equations  written  as  a  first  order  system.  The  scheme 
has  desirable  properties  for  shock  calculations.  For  the  potential  flow 
approximation,  we  show  that  the  entropy  condition  is  valid  for  limit  solutions 
and  that  there  exist  discrete  steady  shocks  which  are  unique  and  sharp. 
Numerical  examples  are  given. 
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SIGNIFICANCE  AND  EXPLANATION 


Hyperbolic  systems  of  conservation  laws  are  often  used  to  describe 
compressible  fluid  flow.  These  hyperbolic  partial  differential  equations  can 
be  approximated  by  finite  difference  schemes  which  in  turn  can  be  coded  for 
computer  calculations  of  practical  problems.  Aerodynamics  is  a  typical  field 
of  application. 

Standard  difference  .schemes  often  run  into  difficulties  when  the  solution 
to  be  approximated  contains  discontinuities  in  the  form  of  shocks  and  contact 
discontinuities.  The  computed  solution  will  typically  either  be  smeared 
i.e.  too  smooth  or  will  contain  unphysical  overshoots  and  wiggles. 

For  a  luge  class  of  scalar  problems  it  has  been  possible  to  design 
difference  schemes  of  upwind  type  which  produce  approximations  of  solutions 
with  shocks  which  are  very  sharp  and  without  overshoots.  In  an  upwind  scheme 
all  differences  are  one  sided  and  the  structure  usually  depends  on  the 
solution  itself. 

This  paper  describes  a  systematic  way  of  deriving  difference  schemes  of 
upwind  type  for  a  class  of  hyperbolic  systems  of  conservation  laws.  Many  of 
the  desirable  properties  which  upwind  schemes  have  for  scalar  problems  can 
thus  be  extended  to  the  physically  much  more  important  case  of  systems. 

This  technique  is  used  to  produce  a  scheme  for  the  potential  flow 
equations.  It  is  proved  that  only  physically  permitted  shocks  will  be 
approximated  in  the  limit  and  that  steady  shocks  are  very  sharp.  Other  cases 
are  investigated  in  computational  examples  which  also  display  the  efficiency 
of  the  scheme. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MRC  and  not  with  the  authors  of  this  report. 
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Many  upwind  difference  schemes  have  very  attractive  properties  when  approximating 
scalar  nonlinear  hyperbolic  conservation  laws.  They  have,  in  particular,  become  the 
standard  technique  for  many  calculations  of  transonic  flow  [1],  [3],  [7],  [10],  [12], 

[16].  There  are,  for  example,  several  versions  of  upwind  schemes  for  the  approximation  of 
the  small  disturbance  equation  of  transonic  flow  (1.1).  A  number  of  these  schemes  have 
solutions  with  sharp  shock  profiles  [1],  [5],  (7),  [16].  The  small  disturbance  equation  is 

(1.1)  2*tx  =  iWx  -1/2  <Y+1)V’x)x  +  V  * 

The  velocity  potential  is  denoted  by  v:x,y,t)  and  K  and  Y  are  positive  constants. 

The  extension  of  these  techniques  to  systems  is  immediate  when  all  the  eigenvalues  of 
the  Jacobian  matrix  of  the  flux  functions  have  the  same  sign.  It  is  the  purpose  of  this 
paper  to  present  a  systematic  technique  for  producing  upwind  difference  schemes  for  the 
more  interesting  case  of  systems  where  the  eigenvalues  of  the  Jacobian  may  have  different 
signs.  We  shall,  as  an  example,  apply  this  technique  to  the  potential  flow  equations  (1.2) 
for  compressible,  inviscid,  isentropic  and  irrotational  flow  [3] .  The  equation  is 

(1.2)  Pt  +  (P^^  +  (P*y>y  “  0  • 

The  density  function  p  is  given  in  terms  of  a  velocity  potential  'P  through  Bernoulli's 
law 
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(1.3)  pY  1  +  (2^  +  *2  +  *2>  ■  H(t)  . 

The  aquation  of  atata  for  tha  prasaure  p  ia 

(1.4)  p  -  ApY 

with  A  abd  y  poaitiva  constants  (1  <  y  <  3  in  our  theorems).  We  can  transform  (1.2), 
(1.3)  into  a  first  order  hyperbolic  system  by  letting  <fix  •  u  ,  -  v  and  then 

differentiating  (1.3)  with  respect  to  x  and  y  resapectively  [2).  The  equality  of  the 
mixed  partials  is  assumed  to  be  valid  throughout  the  flow,  the  system  of  equations  ia 


(1.5) 


Analogous  to  a  standard  procedure  for  scalar  problems  one  might  use  dimensional  splitting 
for  the  solution  of  equations  with  more  than  one  apace  variable,  in  this  paper  we  shall 
consider  the  reduced  one  dimensional  system 
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The  difference  operators  A+  and  A_  denote  the  forward  and  backward  differencing 
respectively  (A±u..  “  ~  u^)).  The  auxiliary  function  f_  and  f+  contains  the 

increasing  and  decreasing  parts  of  f  respectively, 

(1.11)  f.fu)  -  /u  x(«>*,<s)ds 

0 

(1.12)  f  (u)  -  /UU-x(s))f'(s)ds 

0 


n+1  n  At  ...  n. 

Uj  “  Uj  *  S  ^““j1  ‘ 


f  -  f+  +  f_ 

When  f  is  convex  the  definitions  (1.11)  and  (1.12)  simplifies  to 

(1.13)  f+(u)  “  f(max(u,u) ) 

(1.14)  f_(u)  “  f(min(u,u)) 

where  u  is  the  stagnation  point  f‘(u)  ■  0  .  If  f'  has  a  fixed  sign  say  f’  >  0  the 
scheme  reduces  to  the  classical  upwind  scheme 

(1.15) 

When  it  is  used  as  the  basic  Ingredient  in  the  approximation  of  the  small  disturbance 
equation  (1.1)  the  resulting  scheme  will  be  identical  to  the  Cole-Murman  scheme  [16]  away 
from  sonic  and  shock  points.  These  are  the  points  where  f'  changes  sign  and  a  switch 
from  A  f  to  A+f  is  needed.  The  formulation  (1.9)  gives  a  recipe  for  such  a  switch. 

The  approximation  (1.9)  has  several  attractive  properties  in  connection  with  shock 
calculations.  Proofs  and  numerical  examples  are  given  in  [5]  ,  [6] . 

Those  properties  are: 

(a)  The  scheme  is  monotone,  see  [4],  [9],  for  At|f'|  <  Ax  . 

(b)  It  is  in  conservation  form  and  hence  produces  shocks  with  the  correct  location  [14J. 
These  properties  imply  the  following  properties  (c),  ( d)  and  (e)  for  Cauchy 

problems.  (We  thank  P.  Lax  and  M.  Crandall  for  some  helpful  discussions  on  this  matter.) 
In  our  earlier  work  we  proved  the  results  below  for  quadratic  f  even  for  the  mixed 
initial  boundary  value  problem,  and  for  general  scalar  f  for  the  Cauchy  problem. 

(c)  An  entropy  condition  is  valid  for  limit  solutions,  which  rules  out  nonphysical 


shocks 


( d)  The  scheme  is  stable  in  ,  L2  and  LK  . 

(e)  The  approximate  solution  u”  converges  to  u  in  • 

(f)  The  approximation  uses  the  same  number  of  boundary  conditions  as  the  differential 
equation,  i.e.  no  extra  numerical  boundary  conditions  need  to  be  imposed. 

(g)  Shock  solutions  have  stable  sharp  shock  profiles,  [5],  [6].  We  essentially  mean 

that  the  approximation  of  a  steady  Riemann  problem  is  exact  two  points  away  from  the 
shock. 

We  now  try  to  preserve  some  of  these  properties  in  the  approximation  of  systems.  We 
shall  consider  a  hyperbolic  system  of  nonlinear  conservation  laws  in  one  space  dimension 
(1.16)  ut  +  f(u)x  -  0  ,  u  »  R2  *  Rd. 

Systems  in  more  than  one  space  variable  can  be  reduced  to  the  one  dimensional  caBe  by 
dimensional  splitting  or  ADI,  see  [1],  [5]  • 

The  linear  stability  requirement  implies  that  there  only  exists  strictly  upwind 
difference  schemes  if  all  the  eigenvalues  of  the  Jacobian  matrix  3f  have  the  same  sign. 
These  eigenvaluea  are  all  real  since  the  system  (1.16)  is  hyperbolic.  This  stability 
requirement  follows  directly  from  the  domain  of  dependence  and  the  CFL-condition.  The 
upwind  difference  scheme  may,  for  example,  be  of  the  simple  form  (1.15)  if  all  eigenvalues 
of  3f  are  positive. 

We  have  to  clarify  what  we  mean  by  upwind  or  one  sided  difference  schemes  when  3f 
has  eigenvalues  of  different  signs  in  some  region  of  a  solution  space.  Difference  methods 
for  which  the  approximation  of  the  spatial  derivatives  are  non  symmetric  and  may  change 
when  the  signs  of  the  eigenvalues  of  3f  changes  are  often  said  to  be  of  upwind  type,  we 
shall  consider  here  a  special  class  of  such  schemes  for  which  it  is  possible  to  prove  that 
the  properties  (b),  (c)  and  (g)  are  valid  when  approximating  the  system  of  equations 
(1.6).  The  property  (f)  is  valid  in  a  slightly  weaker  form.  Other  upwind  schemes  for 
systems  are  presented  in  [2],  [18],  [19]. 


The  sharp  shock  profile  property  (g)  is  of  computational  importance  and  forces  the 
scheme  to  be  of  a  special  structure.  It  is  otherwise  easy  to  produce  a  linearly  stable 


scheme  of  type  (1.9)  which  approximates  (1.16)  and  is  in  conservation  form.  Choose 
f+  “  ^  (f  +  cu)  and  f_  »  ^  (f  -  cu)  where  c  is  a  constant  such  that  c  >  p(3f).  The 

scheme  is  linearly  stable  for  At/Ax  small  enough.  However,  not  even  for  scalar  problems, 
(d  «  1),  with  steady  shocks  will  the  numerical  solution  have  a  sharp  profile.  This  follows 
from  [11]  since  this  scheme  is  strictly  monotone.  (If  u”+1  =  G*uj+l'uj'uj  i*  then 


3G 


>  0  for  k  =  j+1 , j, j+1 ) . 


The  strictly  upwind  scalar  scheme  (1.9),  (1.11),  (1.12)  is  monotone  but  not  strictly 
monotone  since  at  most  one  of  and  f’  is  non  zero  for  each  u  .  The  natural 

generalization  to  systems  is 


(1.17)  N(3f+(u))  +  N(3f_(u))  <  d 

where  N  is  the  counting  function  for  the  number  of  nonzero  eigenvalues.  This  property  is 
also  basic  for  (f)  to  be  valid.  We  also  need  f  »  f+  +  f_  (modulo  a  constant)  and  the 
matrices  3f+  and  -3f _  must  have  nonnegative  eigenvalues  for  linear  stability. 

When  the  definition  (1.11)  and  (1.12)  of  f+  and  f_  respectively  are  used  in  (1.9) 
the  scheme  can  be  written 


(1.18) 


n+1 


>1 


U j  ”  fx  ^  /n  (1-X<s))f*(s)ds  +  fn  x(s)f'(s)ds 


j-1 


This  is  the  form  of  the  algorithm  that  will  be  used  for  systems.  In  Section  2  we  shall 
present  the  algorithm  for  choosing  the  matrix  x<B>  and  the  path  of  integration  between 
u”  and  u"+1  •  The  choice  of  the  path  of  integration  is  based  on  the  Riemann  invariant 
curves  and  it  is  crucial  for  the  success  of  this  scheme. 

In  Section  3  we  shall  derive  the  explicit  form  of  the  difference  approximation  for  the 
one  dimensional  full  potential  equation  (1.6).  We  shall  also  remark  about  boundary 
conditions  (f)  and  show  that  (1.17)  is  true  for  constant  states. 


Ha  shall  prove  that  the  entropy  condition  is  valid  for  limit  solutions  of  the 
approximation  in  Section  4.  The  global  existence  and  uniqueness  of  discrete,  steady  one 
and  two  shocks  is  also  proved.  These  shocks  are  equal  to  the  analytic  shocks  two  mesh 
points  away  from  the  discontinuity. 

In  Section  S  we  shall  present  results  from  numerical  calculation  with  the  scheme 
derived  in  Section  3. 
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2.  THE  GENERAL  ALGORITHM 

In  this  section  we  shall  present  the  derivation  of  the  upwind  algorithm  for  systems  ?n 
some  generality.  See  also  [17]  for  a  related  discussion  of  this  general  technique. 

Consider  a  strictly  hyperbolic  system  of  nonlinear  conservation  laws  in  one  space 
dimension 

(2.1)  ut  +  f(u)x  -  0  ,  u : R2+  Rd,  t  >  0  ,  —  <  x  <  » 

(2.2)  u(x,0)  =  u(x)  . 

The  flux  vector  f  is  assumed  to  have  continuous  derivatives.  The  difference 
approximation  is  defined  as  follows 


n 


n 


j+1 


(2.3) 


un+1  «  un  -  (  /  (1  -  x(u))3f(u)du  +  /  x(“>  <*£(“)**) 

^  ^  n  n 

u .  u.  , 

3  3-1 


(2.4)  u^  =  u(x^)  . 

The  matrix  x(s)  and  the  paths  of  integration  remain  to  be  specified.  Extending  the 
principles  from  the  scalar  appoximation,  the  matrix  x  should  be  the  natural  projection  on 
Rd  onto  the  subspace  spanned  by  the  eigenvectors  corresponding  to  the  positive 
eigenvalues  of  the  Jacobian  matrix  3f.  These  eigenvectors  are  linearly  independent  since 
the  system  is  strictly  hyperbolic. 

Let  T(u)  be  the  matrix  the  columns  of  which  are  the  eigenvectors  of  3f 


(2.5) 


T  (u)  3f(u)  T(u)  “  A(u) 


f  A,  (u) 


(2.6) 


A(u)  =  diag{A^(u)}  = 


A1  (u) 


Xd(u) 


(2.7) 


*1(u)  <  *2<U>  <  **"  <  ^d(U* 


(2.8) 


X(u)  =T(u)  diag{  I/2  +  V2  sign[A  (u)  )}T  1  ( u) 
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(2.9) 


X(u)3f(u)  “  T(  u)  diag{max( X^(  u) ,  0  )  }T  ^u)  «  (3f(u))  + 


(2.10) 


(I  -  x(u))3f(u)  =  T(u)diag{min(Xk(u),0)}T_1(u)  -  <3f(u)f 


The  choice  of  path  of  integration  significally  affects  properties  of  the  scheme.  The 

path  should  be  specified  in  order  to  simplify  the  computations  and  to  guarantee  good 

characteristics  of  the  solution.  The  definition  will  be  connected  to  classical  techniques 

for  solving  Riemann  problems,  [13]  ,  but  will  be  much  simpler. 

Denote  the  path  connecting  u”_1  to  u"  by  T 3  (and  of  course  r3  +  1  connects  u" 

to  u"+1>.  The  n  dependence  in  T  is  omitted  in  the  notation.  The  curve  is 

decomposed  into  d  subcurves  O'?}?  . 

Jc  K*  1 


(2.11 ) 


d 

u  r3 
k-1  k 


These  subcurves  are  related  to  rarefaction  solutions  and  are  defined  through 

.(k) 

I  It  I 

<  S  <  0 


(2.12) 


r3  ! 

k 


du 

da 


V,kl> 


0  <  s  <  s,  or 
k 


k  =  1 ,  •  •  •  ,d 


where  r^lu)  are  the  right  normalized  eigenvectors  of  3f(u)  corresponding  to  the 
eigenvalues  l^(u).  The  curves  are  connected  by  continuity  conditions 


(2.13) 


( d)  n 

u  ‘  uj-1 

(k-1)  (k),  . 

u  »  u  (s^) 

(1),  ,  n 

1  3 


d,  •  • • ,2  . 


Note  that  the  curve  T3  starts  at  u"  1  with  corresponding  to  the  eigenvector  rd. 


Then  T3  continues  with  r3  etc. 

a- 1 
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The  existence  and  uniqueness  of  a  solution  to  (2.11),  (2.12),  which  is  the  existence 
of  r°  is  guaranteed  if  u"  and  uj-i  are  not  to°  ^ar  aPart*  This  follows  from  the  fact 
that  the  vectors  r^lu),  k  «  1,  •••,<!  are  linearly  independent  and  that  r^  depends 
continuously  on  u  .  In  other  words  TjJ  ,  k  -  1,»»*,d  locally  acts  as  a  coordinate 
system  for  R^. 

An  important  property  with  this  choice  of  path  is  that  the  system  decouples  in  the 
following  sense 


(2.14)  /  din  »  /  x(u(s)  )3f(u(s)  )rk<u(s)  )ds  =  /  max(  A^(  u^(  s) )  ,0 )  rfc(  u(  s) )  ds 

This  follows  from  (2.9)  and  (2.12).  The  scheme  (2.3)  can  hence  be  rewritten 


(2. IS) 


The  eigenvectors  and  the  curves  have  in  many  physically  important  problems  a 

simple  analytic  form.  This  is  e.g.  the  case  for  the  full  potential  equation  and  the  Euler 
equations. 

We  shall  end  this  section  by  showing  that  the  scheme  (2.3)  is  in  conservation  form  and 
is  first  order  accurate. 

uf1  -  u"  “  -  —  (  /  (1  -  +  /  X  (u)3f  (u)du) 

j  j  r3+1  r3 

=  -  "p  (/  <3f(u))  du  +  /  (3f(u))  +  du) 

r^1  rj 

-  -  H  (a_f(u!  +  &+  /  ( 3f  ( u) )+du) 

r3 

=  -  (d+  f(u)  -  A  +  /  (3f(u))'du) 

=  "  Zx  ^2  *0  f(u)  +1/26+  /  1 3f  ( u)  | du) 


-9- 


3.  THE  POTENTIAL  FLOW  ALGORITHM  AND  REMARKS  ABOUT  BOUNDARY  CONDITIONS 


We  shall  now  derive  the  difference  approximation  for  the  one  dimensional  full 


potential  equation 


pt  +  (?*x)x  ■  0 


where  p  Is  a  known  function  of  W  defined  through  Bernoulli's  law 

i  2  c2 

(3.2)  *t  +  '  H(t> 


with  H(t)  given  and 


c2(f>)  “  ^  * 


The  equation  of  state  for  the  pressure  is  given  to 
(3.4)  p  -  A pY  ,  A  >  0  ,  1  <  Y  <  3  . 

In  view  of  the  general  procedure  outlined  in  the  introduction  and  Section  2,  we  shall 
first  treat  the  set  of  equations  as  a  hyperbolic  system  for  the  two  unknowns  P  and 
u  ■  V  .  After  constructing  this  approximation,  we  shall  in  a  future  paper  obtain  a  scalar 
difference  scheme  approximating  (3.1)  and  (3.2). 

By  taking  the  space  derivatives  of  (3.2)  and  using  equality  of  the  mixed  partials 
for  ?  (which  is  assumed  to  hold  even  across  discontinuities),  we  construct  a  first  order 
strictly  hyperbolic  system  of  conservation  laws.  (See  also  (1.6)  in  the  introduction.) 


(  p  )  + 
v  u 


2  2 
£_+  £ — 

2  Y-1  x 


The  2  by  2  Jacobian  matrix 


has  distinct  eigenvalues 


and  right  eigenvectors 


*1,2  “  U  *  ° 


-P/c  p/c 

r,  "  (  ).  r2  “  (  )■ 

1  1 
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Along  the  rarefaction  curves  defined  by 

d  <  P  1  «  f-P/C'i 


(3.9) 


:  :  2-  [  p  }  »  rp/c)  -  r 
1  dt  1  u  ‘  1  1  ;  1 


r  ,  <L  (  P  )  .  (  p/=)  .  r 

2  dt  *-  u  ;  v  1  ;  2 


the  corresponding  Riemann  invariants  are  constant: 


(3.10) 


2  c 

R,  “  u  +  — 7  ■  constant  on  r, 
1  Y-1  1 


2c  _ 

R„  =*  u  -  — 7  “  constant  on  r_ 
2  Y-1  2 


To  find  the  path  of  integration  T3  ,  we  first  find  the  point  of  intersection  In 
the  c,u  plane  of  the  lines 


(3.11  ) 


L  .  u  +  2s_  .  u  +  i 
1  Y-1  j  Y-1 


r  2c  f!i 

L2  !  U  -  FT  ‘  "j-1  ’  re¬ 


calling  this  point  ( Cj  _  i^  #  uj  _  »  **  find 


+  2=1  (u  .  u  , 


(3.12) 


W  2  '  4  '“j  “j-1 


Uj-V/H  <cj  -  Vi* +  V^1 


We 


need  c j  _  >  0  for  the  scheme  to  make  sense .  Thus  our  only  requirement  is 

2 


u.  .  -  u.  <  *7  (c.  +  c.  .)  in  order  that  our  scheme  exist. 

3-1  3  Y-1  3  3-1 

In  order  to  construct  the  scheme,  we  need  to  parametrize  and 

always  do  this  using  c  as  the  parameter. 

Along  r2  we  have: 


We  shall 


(3.13) 


dp 

dc 


£ 

2  r  c 


7  C  T )  -  tSt 

Y-1 
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Hence : 


P(c) 


(3.14) 


/  (3f(p(c)  ,u(c) )  )  +  dc  -  ^max(u(c)  +  c,0)(  *  )dc 

rl  ~ 


1-1 


P(c) 


!  3‘^max(uj-i  "  7=7  CJ-1  +  y=T  c'0)(  1  )*»  • 


'1-1 


Similarly,  along  ,  we  have: 


(3.15) 


7-1 


hence 


(3.16) 


-P(c) 


/  (3f(p(c)  ,u(c) )  )+dc  »  -  —r  /  max(u(  c)  -  c,0)  (  ®  Jdc 

A  *  _  * 


"j  Ci-V, 


-P(C) 


-  t  /  3  B*x((u;j  *  ^rr  cj  •  c>'0,(  °  )<*c 


'i-’/i 


In  order  to  carry  out  the  integration  in  (3.14)  and  (3.16),  we  first  use  the 
indefinite  integral  results: 

P(c) 


(3.17) 


and 


-2-  /  (U  -  -i_  c  ♦  HI  C)f  C  )dc 

Y-1  1  '  j-1  7-1  j-1  7-1  H  1  J 

2P<£l  c 

■«v--^Tci-,,(r 

7- 1  C  2  ° 

T  (7-1) 


(3.18) 


-P(c) 


7=7  /  (uj  +  7=7  cj  "  iM  c>(  i  J 


dc 


2 

P(c) 

•f  l 

7=1  P(C)  C' 

“  (u,  +  — -  c. ) 

1  Y-1  l' 

-2 

7+1  2 

Y-1  ° 

2  C 

(Y-1) 

t 
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L 


We  may  use  these  expressions  in  order  to  evaluate  the  inte9rals  in  (3.14)  and  (3.16), 
arriving  ats 


,  ,  P ^  lx  —  P  .  * 

/  (3f(p(c)  ,u(c) ) )  dc  -  (uj_1  -  cj_1)(  /2  3 


2  ~ 


(3.19) 


Y=T(Cj- V  cj-i’ 


) 


2  ~ 


TM  (P)-%C)-V2-  °3-1C3-1) 

+  (  2  ) 
-1+1-  (e  ,  .  ~2 

l-’4 


where  p  «  p  ( c)  ,  and 


(3.20)  (a)  if  min(u;j_1  +  c^_t/  “j_1/2+  cj- *  0  then  cj-1^“  c j-  1^'  cj-l  "  Cj-1 


(b)  if  uj.,  +  c^  >  0  >  u^.  1A+  Ci.  1 


j-  ^  T  Cj-V2 


then  VV  *j-V "  ni  (Vi  ‘  wca-i*  and  Vi  “  cj-r 


(c)  if  u_j_.,  +  Cj,,  <  0  <  Uj_1^  +  cj_  V2  then 


W  cj-^  and  cj-i 


(^)(uj-i  *  pr  cj-i} 


(d)  if  max(uj.1  ♦  Cj.,,  Uj_  1^+  Cj_  l/2>  <  0  then  Cj_  cj_1  -  Z ^ 


(and  in  this  case  the  quantity  in  (3.19)  vanishes). 


We  have  defined  c 


.  i,  so  that  on  it  is  true  that  u.  ,  »  -c.  i.e. 

j-  Vi  2  3-  Vi  j-  v2 


it  is  a 


sonic  point 


M«  AlSO  hAVSI 


/  Of  (p(c)  ,u(c) )  )*do  -  -(u  ♦  c  ; 

4  J  T“*  J 


(3.21) 


rJ 

i 


PJ  '  P3-  ^ 

"  7h(ci  ■  CJ- 


2  - 


■  Tf-1  <P3Cf  Pj-Vfj-W 

HJ  ,;2  -2 

2  C1  ci-  V» 
(X-1)3  3  3  * 


whsra  p  -  p(c)  and 


(3.22)  (a)  it  ain(Uj  -  Cj,  “j.  ”  cj- ^  *  0  **>•"  CJ“  Vi'  °3  *  °3 


(b)  it  ttj  “  cj  *  0  >  “j-1/*,-  CJ- ^ 


th,n  SJ-V  <’V£ei>'  V  c3 

(c)  if  Uj  -  Cj  <  0  <  Uj_  V  than  ^  -  Zy  ^  ^  Cj_  ^ 

N  —  — 

(d)  if  H«(Uj  -  Cj,  Uj_i^-  cj- 1^  }  <  0  than  c j  ■  CJ- 

(and  in  this  cans  ths  quantity  In  (3.21)  vanlshss). 

Wa  havs  dafinad  c  ao  that  on  rj  it  is  trus  that  u  1  »  c  l.a.  it  is  a 

}~  ri  i  3“  fi  j~  fi 

sonic  point. 

Using  dividad  diffaranca  oparatocs  wa  nay  now  writs  out  tha  axplicit  ona  stap  upwind 
dlffaranca  approximation  to  (3.8)i 


(3.23) 


X 

n  n 
°i  u3 

-  -D+ 

n 

,  n.2  ,  n,2 

lU3 

(Uj)  (Cj) 

„  2  +  y-'  . 

.  _x_  ,  n  n  n  n  . 

+  D+A+(Pj,UJ,Pj.1,Uj.1)  - 


-  -  JZ  /  Ot(wn))+«Jw  -  ~  J  (3f  (wn)  )”dw 


Ax 


r3>1 
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V ' 


In  order  to  invert  the  peir  of  non-linear  aquations  for  (  ^  ) ,  we  suggest  using 


P 

Newton's  aethod  with  the  initial  guess  (  ^  )  ,  the  solution  at  the  previous  tine  step. 

U" 

It  is  suggested  by  several  authors,  e.g.  [1],  [20],  that  only  one  iteration  suffices  in 
cases  such  as  this. 

Thus  the  full  expression  ia> 


(3.28) 


n+1  n 

(  l,  .)■'((  i  )*3ir 

p  ,U  Uj  J  J 


Of  course  in  the  regions  where  the  Jacobian  matrix  has  eigenvalues  which  are  both 
positive  or  both  negative  the  matrix  dTp  u  is  aither  upper  or  lower  triengular.  In 
particular,  if 

than 
(3.29) 

and  if 


them 

(3.30) 


Uj+r  >  Cj+r  ,  r  -  0,  *  1 


dT  (  ^  1  •  -A  ( 

£.u  v  u^  1  -  v 


*  -uJPi 


) 


vj  +  rr 


-J+r  <  fj+r  *  r  -  0,  *  1, 

f  p,  .  f  £iui *  Vi  . 

dTp,u  ^  u  .  ^  ^  2  ^ 


Vj  *  £j  PJ 


In  general,  a  complicated,  but  fairly  straightforward  calculation  gives  us  the 


linearised  operator  at  a  state  £^,u^ 


''  *  *- 


(3.31) 


dT  (  M  “  “D*( 

p,u  1  '  +v 


£jUj  +  Vl 


)  +  D*[max(u^  -  cjt0)(  _c  ) 


=* +  ir p* 


Si  P3 


t“j-v2- 3_t  ^1/2^pj  +  l^;pj.i  +  uj-u3-i1 


+  (  -j  -J-  'h.  )(U  +  rl  pj+nuuciu  ,  +  c  '/2j  V2  Cr^  PS  p,., 

.  _4  (=  1()  3  £J  j  j  4  3  *  ’  P-i  3  £i-’  3 


r-T  -3  Cj-V2 


1  £j~  V2‘  £j-i  -j-1 

-  c  *  (  J_(~  -  ><V.  -  £- 

-J-L  Y-i(Sj-V2  £j-D  j  1 

In  the  subsonic  case  whan  the  eigenvalues  of  the  Jacobian  matrix  3f  are  of  both 

signs,  the  operator  dT  Is  a  perturbation  of  a  simple  form,  in  particular  if  -  c.,.  _ 
P  »U  “JT 1 

jij+r  <  £,j+r«  r  "  0,  *1  ,  then  we  define 


(3.32) 


dT  (  Pj  )  -  -  I  B<J)<P.u>(  P)+r  ) 

o.uv  u  y  r  •  —  v  u  7 

-  -  j  r— 1  r  Uj+r 


prove  the  followingi 


B_j'(p,u)  -  B_t  (£^_1 ,Uj_1  )  +  0(|pj  -  Pj_1 I  ♦  lUj  -  U  I 


8  ^  (pell)  *  B  (p  „  #u  )  +  0(|p  —  p  I  +  |u  -  u  I) 

1  --  1  - j+1  '-j+1  '  j+1  1  'j+1  j 


(a)  B_1  and  B^  ara  both  of  rank  ona. 


(b)  Tha  aiganvaluaa  of  B_1  ara  0  and  _Uj_1  +  Cj_1  >  0  ,  with  corresponding 

alganvactora  (  ). 

1  1 

(c)  Tha  aiganvaluaa  of  B^  are  _u^+1  -  _c^+1  <  0  ,  and  0  and  tha  corresponding 

ij+1  £jii 

c  c 

eigenvectors  are  (-  )  ,  (  ”^+1  )  . 


As  a  consequence  of  this 


we  have  the  followings 


Rssark  (3.1) 


At  a  boundary  point  near  which  the  flow  is  smooth,  tha  affect  of  nonphysical  and 
purely  metrical  boundary  conditions  la  small  (in  fact  saro  In  tha  aupersonlc  case).  This 
is  true  because  of  the  previous  Lasts  since 

X(P.u)3f (£<u)  -  B_? (p,u) 

(1  -  x<£«u>  )3f  (p,u)  -  B^p.u) 

(Bee  <3. 8)  for  the  definition  of  x<£>£))  • 

Proof  of  Least  (3.1) 


We  use  (3.31)  In  this  region  to  write 


(3.33) 


>  V  Si-, 


1  c 

<“3-,  *  Sj-l’l  c  ^  ^  ^ 

3  3  -I-1  £J-1  2  .  . 


7^(-cJ- V  fj-l) 


Si-  u 


-1-  ’A 

♦  '‘Sj-’A*  Sj-lJ  ♦  ( 


£3-V£3-l 

J 

r^T  <£j-  V2"  £J-i’ 


UPt"'-.  M 


So 


(3.34) 


j-i  _  -j-i 


Kl-1  =1-1  -1-1 

i,  (P,  ,)(  ;  »  P  [(  -  («!_,  + 

-1  -j-1  -j  1  Uj_i  j-1  1  £j.1 


*  Uj_1l(  )<’/2<“j.1  +  fj-i))] 


and  hence,  finally: 


(3.35)  B_1(Pj_1  ,£._,)  -V2  (Uj.,  +  5 j-1 * 


-j-1 


g-j-1 

-j-1 

1 


A  simple  calculation  shows  us  that  u^,)  has  one  nontrivial  eigenvalue 

ej-1/£j-1 

equal  to  u^,  +  ,  and  corresponding  eigenvectors  equal  to  [  ^  J  .  The  zero 

"£  j-1 /-j-1 

eigenvalue  has  corresponding  eigenvector  (  ^  )  • 


An  analogous  proof  works  for 
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4.  THE  ENTROPY  CONDITION  AND  STEADY  DISCRETE  SHOCKS 


Our  first  theorem  of  this  section  is  thus: 


Theorem  (4.1) 

P^t) 

u  j  ( t) 

to  w(x,t)  with  inf  p^(x,t)  >  >  0  .  In  addition,  suppose  the  quantity 

lim  sup  ( | A  p  . ( t) |  +  |A  u.(t)))  is  sufficiently  small  (but  positive).  Then  w 
Ax*0  j,t  3  3 

satisfies  the  entropy  inequality  (4.5). 

Next,  we  present  analytic  evidence  that  our  scheme  gives  excellent  shock  resolution  in 
the  steady  case.  (Numerical  computations  are  presented  in  the  next  section  both  for  the 
steady  and  unsteady  case . ) . 


j  solves  (4.6)  and  converges  boundedly  a.e.  as  Ax  ♦  0 


Suppose  w^(t)  «  ( 


L  R 

p  p 

Let  (  )  ,  (  )  be  the  states  on  the  left  and  right  (of  x  =  0  say)  for  a 

L  R 

U  u 

steady  shock  solution  of  (3.5).  This  means  that  both  the  jump  conditions 

L  L  R  R 
p  u  »  p  u 


(4.7) 


,  L,  2  ,  L  2  ,  R  2  ,  R,2 

(u  )  (c  )  ^  (u  )  (C  ) 

2  Y-1  2  Y-1 


and  either  the  conditions  for  a  one  shock: 

(4.8)  uL  >  c**  and  -cR  <  uR  <  cR 
or  two  shock: 

(4.9)  -c1  <  uL  <  cL  and  uR  <  -cR 
are  valid. 

Then  we  seek  solutions  of  (3.23),  (3.25)  or  (3.26)  which  are  time  (or  n)  independent 


and  approach  (  )  at  and  (  )  at  +•  .  These  are  called  steady  discrete 

L  R 

.  .  u  u 

shocks . 

We  have 

Theorem  (4.2) 

(i)  Existence:  A  class  of  steady  discrete  shocks  exists  globally  and  each  member  is  of 


-2  1- 


Here  the  solution  of  the  system  of  differential  equations  is  to  be  taken  for 


0  <  a  <  a  ,  with 


(  )  -  ( f  ) 


V’"’’ 


,  “  Y-1  ,  R  2  R. 

and  c  -  (u  +  —  c  )  . 

(b)  2  shock: 


L  ,  2_  ,  ,  .  L. 

u  “  u  +  7T7  *c-t  ta)  -  c  ) 

Jfi  ‘  '  Jn 


-  -  L  L 
-p  u  +  p  c  +  p  u 

■'0  30 

V1  - 


(a) 


(cj  +1<«)  -  «j  +1<“>) 
2  0  0 


V  ^  eV,<“> 


cs  +1(“)  -  ~ r  <C  (a)  +  U  <C0>( 

Jn  '  1  Jn 


W'"11 

W'"" 


- 1) 


with 


c  (C) 


V’(0) 


)  -  (  ) . 
c 


Here  the  solution  of  the  differential  equation  is  to  be  taken  for  0  <  a  <  a 


\»u*‘ 


)  - 1  ) . 

R 


.  *  r  111  i  ,  L  2  l, 

and  c  -  -l  ^77  j  (u  -  —  c  ) 


with 
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(ii)  Uniqueness: 


We  consider  only  steady  discrete  shocks  having  a  certain  weak  monotonicity  property. 

J  j/2  +  °j/2 

all  integers  j  ,  and  (ii)  if  u^  <  c_j  then  u_j+1  <  Cj+1* 

(b)  Admissibile  discrete  steady  two  shocks  are  such  that  (i)  u ,  -  c.  <0  for  all 

j/2  ]/2 


and  (ii)  if  u^  <  -c^  ,  then  u^+1  <  -=j+1  • 


intergers  j 

Under  these  assumptions,  all  discrete  steady  shocks  are  of  the  form  given  in  (i). 

We  note  that  without  any  of  these  restrictions,  discrete  steady  shocks  must  be 
eventually  constant.  This  means  there  exists  integers  a  <  b  such  that  E  /  of 
j  <  a  ,  and  w^  =  w  of  j  >  b,  for  some  a,b  .  This  is  the  content  of  Corollary  (4.1) 
below. 


We  also  believe  that  our  uniqueness  result  is  valid  under  weaker  hypotheses  -  see  the 
scalar  result  in  [6] ,  [7]  . 

Proof  of  Theorem  (4.1) 

In  its  broad  outline,  this  proof  will  follow  the  entropy  inequality  proofs  in  [6], 

17).  There  are  various  differences  due  to  the  vector  valued  nature  of  this  problem. 

Let  Wj(t)  satisfy  (4.6)  and  let  ^(x,t)  >0  be  in  cJ(R+,R)  .  Multiply  (4.6) 

by  ^(x.,t)  VT(w,t)  Ax  ,  sum  and  integrate,  and  add  T  /  *>(x  .,t)DXG(w  ,(t)  )dt  to  both 

3  w  j  3  +  3 

sides.  We  then  have 

(4.11)  -/  dt  l  Ax[^t(xj,t)V(w;j(t))  +  (D^x)v>(xj,t)  JGtw^t))] 


J  dt  £  Ax*(x  ,t)[  /  VT  (w)[(3f(w))+  +  (9f(w)) 
j  3  r 

j+i 


dw 


-  /  V^(w.)((8f(w))*)dw  -  /  vr(w)((3f(w))+)dw)  . 

p  w  3  r  w  J 

3+1  3 


As  Ax  ♦  0  the  left  side  converges  to 

-  / /  (  V^V(w)  +  V^Gfw) )  dxdt  , 


(4.12) 


by  the  Lebeague  dominated  convergence  theorem.  We  must  merely  prove  that  the  right  3ide  of 
(4.11)  is  nonpositive  as  Ax  +  0  . 

We  rearrange  the  terms  above,  drop  the  t  dependence,  and  arrive  at 

(4.13)  Right  side  of  (4.11)  =  / dt  i  A x*>(x.)  /  [VT(w)  -  VT(w  ) ) ( 3f ( w) ~) dw 

j  3  r  ^  ^  j 

3  j+i 

+  f  dt  I  Ax^(x. )  /  (VT(w)  -  VT(w  )](3f(w))+dw 
‘  “  3  i,  w  w  j+1 

3  j+1 

+  /  dt  l  Ax(D*  v>(x  ))  /  V^w  >(3f(w))  +  dw  -  [I]  +  till  +  [HI]  • 
j  3  w  3 


It  is  clear  that: 


|  [III] |  <  K  /  dt  l  Ax  |A+w  |*x 


where  K  depends  on  •f  .  Hence  [III]  ♦  0  by  the  Lebesque  dominated  convergence  theorem. 
Next  we  consider  the  integral  along  I*  ^  in  [I] 


(4. IS) 


/  <VT(w)  -  VT(w.)  )3f(w)"dw 

rj+1  ”  "  j 


/  (VT(w)  -  VT(w.))(3f(w))+dw  +  /  (VT(w)  -  VT(w  n(3f(w))'dw  . 

J.j+1  W  “  3  J.j+1  W  “  3 

2  1 


(4.16) 


/  (VT(w)  -  VT(w.))(3f(w))'dw 

rj+i  W  W  j 

2 


=  -~r  I  j+1/2dc  min(u - c  +  c,0)[-^- 

Y-1  ‘  j  Y-1  j  Y-1  Y-2 


2  c .  ,  . 

H - JL  )£!£>.  +  _2_  (c  .  c  ) 

P ( c )  '  c  Y-1  j' 


l-h)2  JCj+V2  do  min(u  --^c.+^c,0)  [JC  l  ~fr  —  +  l)dv] 
lY-1  '  j  Y-1  3  Y-1  1  P(v)  c  ‘  J 
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r— ~ 


Next  we  have  ■ 
(4.17) 


/  (V*(w)  -  V*(v  ))(3t(w))~dw 


r3+’ 


’  I  3+1min<Uj+1  *  ^7  °j+r  £l  C'0)[^2( 
CJ+V2 


2  c.‘  ,  , 

JL  ifilsl 

P ( c)  ’  c 


"  uj  +  y=T  <cj+i  "  c>,dc 


^  /  3  -in(uj+1  ♦  ^  cj+1  -  c(0)[^  - 

C4+  U  ^  /2 


i*  V4 


2  _2 

c„  1,  -C, 


(c .  Cj+  ^  +  _i_  (  _±l26  -i  )  £i£i-  £  (Cj+  ^  Cj)]dc 


v _  p(c) 


-  !  3+1  dc  maxlVi  °3+i  lp(v>  c 

CJ+  Vfe  c  j+  v2 


+  1)dv 


+  cj’  ] 


where  K  and  K  below  are  uniformly  continuous  functions. 
Similar  analysis  gives  us 


(4.18)  /  (VT(w)  -  VT  (w  ))(3f(w))+dw 

rj*1  w  3 


-  -(£■  )aJ  * *  ~X‘»J  -  *  •,  ♦  xti  c,o,[  J3+  etal  +  , ,dv 


+  *‘C'cj+,'cj*V5><cj+i  "  cj+i/] 


and 
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(4.19) 


/  <V*(w)  -  V^(wj+1))(3f(w)>+dw 


-  ^)2  ”“«vi +  y=t  °j+i  -  ?r  c'0)  f  j+1(^rr  421  ♦  ,)dv 


c3+  v2 


He  add  the  last  four  equations  and  arrive  at: 


(4.20)  /  [vT(w)  -  V*(w  )](3f(w))"dw  +  /  [vT( w)  -  VT(W  )](3f(w))+dw  < 

•I..  w  "  J  J.,  H  W  3  +  1 


.3+1 


3+1 


<  - 


.(1)  .  .(1) 


(2)  ^  .(2) 


- -  [-min(u  +  c  ,0)  +  max(u  +  c  ,0)](c  -  c  ) 

( Y-1 )*  J  J  J  3  J+  *2  J 


(2)  12)  2 

+  a  max(Uj  +  Cj  ,0) |cj+  ^ -  c J  |cj+,  -  cj+  ^ 


(Y-1) 


^~2  [-min(u*3)  -  c*3,,0)  +  max(u*4)  -  c*4)  ,0 ))  (c.,+1  -  c^+  1a)2 


3  3  . . 3+1  3+ 


•  *  "  cj4>,0>  |cj+  i/Cj“  c3'2'cj+,  ■  C3+1^  • 

Here  a  Is  a  universal  positive  constant  and  the  u^1^  <  are  evaluated 

somewhere  on  the  relevant  above  paths  of  integration. 


We  first  note  a  simple  fact: 


-mintu*1  ^  +  c^,0)  +  M«(u^  +  c*2*,0) 


(4.21) 


-min(u^3^  -  c^3*,0)  +  maxtu^4^  -  cj4^,0)  >  &  >  0 
for  6  a  universal  constant. 

The  right  side  of  (4.20)  is  obviously  nonpositive  for  |A+w^|  <  6  *  except  perhaps  ir 
the  following  cases:  either 


(a)  inf  |Uj  -  y!?  Cj  ♦  £1  Cl  <  =up  |ujtl+^j+)-gc|| 


°3+1  '  C3+  V2' 


c  e  r_ 
or  2 


3+1 


c  €  r 


3+1 
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* 


1 


(b) 


infj+1lttj+i  c|  4  aup4Juj  ’  y^c3+£tc||c3+V  cj' 


c  e  r 


c  £  r 


j+i 


1  2 
However  when  (a)  happens,  it  is  easy  to  dominate  (4 .20)  for  e  sufficiently  small, 

(4.22)  /  |VT(w)  -  VT<w.) |(3f(w)>"dw  +  /  |VT(w>  -  V^(w  ) | ( 3f ( w) ) +dw 

w  w  3  '  ...  w  w  j+1 


rj+1 


j+1 


4  -bi|c3+V  °3 |3  "  b2|o3+i  "  W/  +  b3|c3+V2‘  °3|2|c3+i  ‘  °3+  V/ 


+  b4|c3+ V  C3I|C3+1  ‘  °3+W 


for  the  bj  universal  positive  constants. 

It  is  easy  to  show  that  this  expression  is  nonpositive  for  sufficiently  small  e  . 
A  similar  argument  follows  when  (b)  occurs.  Thus  the  expressions  [I]  +  [II]  is 
nonpositive  and  we  are  finished. 

Proof  of  Theorem  (4.2) 

He  begin  by  demonstrating  that  the  functions  in  (4.10)  are  Indeed  steady  discrete 
shock  solutions  of  (3.23),  (3.25),  or  (3.26).  This  means  that  the  sequence  Wj  * 

P3 

(  J  solves  the  following)  for  each  3 

Uj 

(4.23)  /  Of  (w)  )*dw  +  /  Of(w))”dw  -  0  . 

rj  rj+1 


This  is  trivially  valid  for  3  <  3Q  ”  2  ,  3>3fl+3. 


we  must  merely  verify  that 
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(4.24)  (a) 


(b> 


(c) 


<d) 


J.  (3f(w))+dw 

r  0 

/  +  (3f(w))+dw 

r  0 


/  (3f(w))"dw  - 

r  0 

+  /  ,  (3f(w))  dw  ■ 

r  0 

+  Jj  +2Of(w))'dw  - 

r  0 

/  .  +,<  3f(w) )+dw  ■  0 
'n  £ 


0 


0 


0 


We  need  only  verify  any  three  of  these,  the  fourth  will  then  be  valid  automatically. 
This  follows  by  merely  summing  all  four,  and  using  the  shock  jump  relation  f(uL)  •  f(uR). 
We  first  consider  the  discrete  one  shock  case. 

In  order  that  (4.241(a)  be  satisfied,  we  need  u  >  -c  on  the  line  segment 
connecting  (cL,  uL)  to  (c.  i,,  u,  ..)•  (We  shall  always  require  that  c,  i,  >  0 

i0  ~  '2  j0  -  h  30  +  >2 

for  v  -  -1,0, 1,2,3  so  that  cavitation  does  not  take  place.)  we  also  require  that  u  >  c 

on  the  line  segment  connecting  (c.  j.  ,u.  j.)  to  c  ,u.  ).  This  makes  the  condition 

' q“  '2  ' q~  '2  Jq  'o 

above  redundant.  Thus  for  (a)  we  need  only  u,  >  c.  >  0  ,  u,  ,  >  c  , ,  >  0  . 

10  J0  J0-  '2  J0"  '2 


For  (d)  to  be  valid,  we  first  require  that  (c  ,  u  >  -  (c  +  ,u  +  )  . 

•  Jo  Jo  Jo  '  ' 

This  means  that  (c  ,  u  +1 )  is  connected  to  (cR,uR)  via  a  one  wave  -  i.e. 

■'O  -*0 

2  R  2  R 

u  +  +  —  c .  -  u  +  —  c  .  Finally,  for  (d),  we  require  u  <  c  >  u 

■'n  '  ■'n 


We  now  must  verify  equation  (c).  Since  (c  ^  ,  u  +  )  is  connected  to  (cR,uK)  by 
a  subsonic  one  wave,  we  have: 
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(4.25) 


/  +1Of(«))"<Jw  -  f(wB)  -  f(w  ) 

j.30  30 


p,luR  ■  V1  V1 


R.2 


(c  ) 

t-i 


*>2  V1  V1 


Jo 

T-1 


Next  we  require  that 

(4.26) 

where 


u .  .  >  c .  ,  ■).  >  0  ,  which  implies,  among  other  things! 

30  ^  V  '2 

1 4  ,  (3£(w)  )+dw  -  f(w  1z)  -  f(w  ) 

3o  1  3n 


vv  (W  V-4* 


(  Sr  ,<uj0+i +  tS  cj0+i)(1'1> 


T-1  ,  R  .  2  R,  .. 

“  7m  (u  +  To  0  ’  (,'1) 


(ass  aquation  (3.2T) (b) .) . 


Thus  to  verify  that  wa  indeed  have  a  steady  discrete  shock,  we  need  to  show  that 


(4.27) 

ori 


f(w")  +  f(w„  .,)  -  f(w.  )  -  f(w.  _.,)  -  0 

30+  ^  3o  30  1 


(4.28) 

V 

\  V  ‘  Pi0\  * 

R  R 

u  -  Pj^ 

V1 

0 

-2 

V 

’c  v  v2  \ 

2 

cj«  ,  R»2 

Jo.  +  to 

,  (cV  V_ 

2 

V1 

* 

2 

+  T-1  2 

T-1  2 

T“1 - 2 

T-1 

V 

subject  to  the  restrictions  that 


(4.29)  (a) 


u .  >  o  t,  >  0 

V  V2  3n"  V2 


u  >  c  >  0 

]0  30 


V  ’4  *  V  v2  '  0 


—  C  (  u  <  c  >0 

j +i  V+i  j  +i 

0  0  0 


u  +  -2-  C  -  uR  +  -2-  CR  . 

V1  Y-’  V1  Y_1 


Next  we  notice  that  the  jump  conditions  (4.7)  are  equivalent  toi 


(4.30)  (a) 


(PRuV  ,  (cV  (uR)2  ,  (cR)2 

2<pV  7-1  2  7-1  ' 


r  «  R 
L  p  u 

"  mtF 


The  derivative  of  the  function  on  the  left  in  (4*30) (a)  with  respect  to  c  is 
(PLCL)  2((pLcL)2  -  (pRuR)2).  Thus,  a  single  minimum  occurs  when  uL  -  cL  - 

7-1  _2 

R  7+1  R  Y+1 

(u  ) 1  (c  )  .  For  this  value,  the  quantity  on  the  left  in  (4.30)  is 


—  2(i=i) 

+  1 _  .  R.7+1  ,  R.  '7+V  „  1  R  2  i.  .  R  2 

pyy  (c  )  (u  )  <  —  (c  )  +  V2  (u  )  , 


by  Holder's  inequality#  Thus,  the  jump  conditions  have  one  solution  (c^u")  for  each 
(cR,  uR),  and  for  0  <  uR  <  cR,  we  have  uL  >  c1,  >  0  .  This  solution  is  characterized  by 
the  root  of  (4.301(a)  satisfying 

2  7-1 

„  „  L  ,  ,  R.7+1,  R.7+1 
0  <  c  <  ( c  )  (u  > 


This  analysis  is  needed  in  what  follows. 
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Solving  the  first  aquation  In  (4.28)  for  u  glvaa 

30 


(4.31  ) 


R  2  .  R  .  .  R  R  ®  . 

pj  *i(u  +  y^T  <c  *  cj  ♦i”  +  p  u  +  p  c 

Jo  Jo 


in 


(dropping  tha  j0  *  Vj  subscript,  both  abovs  and  In  what  follows). 

The  second  equation  In  (4.31)  than  becomes 

R  .  2  ,  R  ,1  RR  ---,2 


r  ,R  i  ,  R  ,1  RR  “  • 

lpj  +i(u  (c  ‘  cj  *i 1 J  “  p  u  -  p  c  J 

•'n  •'rt 


o  » 


2p* 


(4.32) 


.  J°  +  _JLU _  (;,2  +  liiii  +  JiLii  .  Vl  . ,  (UR  +  JL.  (c*  .  c 

Y-1  2(Y-1)  ‘  2  Y-1  Y-1  ^  Y-1  j  +1 


)>‘ 


thus  defining  the  function  g(c,  ,  c.  )  . 

30  30  1 

One  solution  to  this  aquation  Is  (c  ,  c  )  “  (c,c)  ,  anothsr  is  (cL,c). 

J()  ^0  1 

shall  show  that  there  exists  a  one  paraaieter  family  of  solutions  (c  (a),c  (a)), 

R  L  _  J0  30 

both  components  increasing  monotonically  from  (c,  c  )  to  (c  ,  c) . 

He  shall  then  demonstrate  that  (4.29) (a)  -(d)  is  valid  for  these  solutions, 
immediate  from  the  construction) . 


K  straightforward  calculation,  using  (4.31)  and  (4.32),  gives  asi 

3  2 


(4.33)  (a) 


3c 


in 


(b) 


u  p 

3c  "  pr  (uj  +i  •  cj  +  p  J  °'  ) 

V1  T  jo  3o  Vv1 


we  consider  the  system  of  ordinary  differential  equations 


We 

with 

((e)  is 
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(4.34)  (a) 


dc  ^  (a) 

"  V1  <Cj° '  Cj°+1> 


dc .  .(a) 


with  initial  conditiona 


c  (0) 
30 


)  “  (  ) 


C.  ^  (0)  R 

V1  C 

to  b«  solved  for  a  >  0  . 

Since  the  initial  data  solves  (4.32),  it  is  clear  that  the  solution  of  (4.34)  will 

also  solve  (4.32).  We  shall  now  show  that  g  >  0,  g  40,  for 

c  c 

h  V1 

m  L  R  m 

c  «  c  (a)  <  c  ,  c  4c  (a)  4  c  ,  with  g  »  0  only  at  the  left  endpoint 
J°  j°  % 


:  “  c  ,  g  “0  only  at  the  right  endpoint  c  •  c  . 

30  Cj  +1  30 

-*  n 


The  statement  concerning  g  is  simple  to  verify  since  by  (4.321(e), 

V  1 


V1  ■  V’  V1  “ c '  and  %  +1 < 0  at  a  ■ 0 


In  order  to  prove  the  statement  concerning  g  ,  we  shall  use  p  u  >  p  c 

CjQ  30  :0  J0 

for  all  the  relevant  values  of  a  except  a  -  0  . 

This  is  equivalent  to  showing:  for  these  values  of  a: 


(4.35)  p  (u  -  c,  )  ”  p  u  +pc-  (p.  c  +  p  u  )  >  0 

Jo  Jo  ]o  3o  V’  V’ 
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Equality  hold*  for  a  »  0 


Using  (4.31)  -  (4.34).  wa  differentiate  the  left  aide  of  (4.3S)  and  obtain: 

P- 

£_p(u  .  c  )  .  .  (  lil  )p  c  --2-e 

do  jQ  jQ  5„  1  Y-1  JP"  « 


V1 


0  Jo 


c,  .  j„+i  j.+r 

“  .+1  0  Jo 


J0  Jn 


2  r  u  2  ,2  2  ^  J0  .  Y+1  , .  .  J0  J0  m 

^ v  ‘  V’![^  i v  -i.’v;  ^  \1’  ^v" 


U  D 

3n  Pjn+1 


Thus  we  have  the  inequality: 


(4.36) 


—  p  (u  -  c  )  +  X  <a)p  (u  -  c  )  <  -K  (a) 
do  Jg  Jg  Jg  1  3g  J  g  jg  2 


with  Kj  >  0  as  long  as  u  <  c  1  . 

20  30 


This  mans  that  u  (a!  >  c  (o)  >  0  as  long  as  u.  (a)  <  c.  ..(a)  . 

Jq  Jg  Jq  ^  ^0  * 

Thus  the  solution  of  (4.34)  sxists  inside  cL  >  c  >  c  ,  c  >  c  >  cR  .  We  wish  to 

L  ■  ■ 

show  it  escapes  at  (c  ,  c)  .  Whan  the  solution  passes  through  (c,  (a  ),  c  )  ,  as  it 

30 

■uet  for  som  a  ,  equations  (4.31),  (4.32)  are  valid  with  u  >  c  , 

m 

U  ■  C  ■  C  e 

V’  V1 

Thsss  squat  ions  then  bscoa#  the  jisnp  conditions,  (4e30),  already  analyzed,  with 

L  L 

(p  ,  u  )  taking  th  place  of  (p  ,u  )  .  Thus  the  unique  solution  must  be 

c  (o^  -  CL  . 

J0 

We  saist  now  mrely  verify  inequalities  (4. 29)  (a)  and  ( c )  .  ((b)  and  (d)  are 
issMdlate  tram  the  construction  1. 

First,  we  note  froai  (3.12),  that 


L 

c  ♦  c 


(4.37) 


:  .  -  - ^  (u  -  uL>  . 

V  4  2  4  10 
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Thus,  using  (4.31)  -  (4.34),  we  have 


dc^  du,  dCj 

_ fo 

da  4  da  da 


dc  du,  dc,  u 

(4.38)  ~  c,  -  V2  — ^  +  IZl  1^-12  .  JO  j  _i 

da  j  -  ^  da  4  da  *  da  v  c  ' 


dCjo+1  V1-  Uj«+1 


da  p . 


-(i-r-2-) 


V1 


i.  ;%■  v ,  W1  p^+i 

~(v  — 


(i 


<c  +  U  )  ) 


PJ0CJ0+1  PioCjo+1  3°  j° 


■  1-'  <Uin+1  "  V1<( 


c  ~  U  CD 

j°—  j0)(i.A  V1 


)n 


\  V1 


)  • 


This  quantity  is  easily  shown  to  be  nonnegative  if  1  <  y  <  3  .  Hence  the  minimum  of 
c.  ..  occurs  at  u.  “  c.  “  c  . 

v  4  10  i0 

We  substitute  that  into  (4.37)  and  we  must  show 


(4.39) 


cL  +  cR  .  r-1  .  R  L,  .  . 

- ^ -  ♦  J-r-  (u  -  u  )  >  0 

2  4 


This  means  we  must  verify: 


(4.40) 


2cL  .  2cR  .  R  .  L  /  R,2  .  2(cV  2U? 

-  +  ——  +  u  >  u  »/(u)  ♦  — -  - - — 


V 


7-1  7-1 


7-1  7-1 


(using  the  second  jump  condition  in  (4.7)).  This  is  trivially  valid  if 
2 


7-1 


>  1  ,  or  a  <  3  . 


(4.41 ) 


Next  we  note,  again  using  (3.12),  that 

1 


U  V  V2  '  7-1  (Cj0  -  CL>  ♦  2 


0  L 
+  U 


and 


(4.42  ) 


dc  .  du . 

d_  1  _Jo  ,  ,30 

da  Ujjj-  Vj*  7-1  da  ^  da 


dC 


V  V2 


r-1  da 
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Thus  i 


(4-43)  'c- 


(Y-1,2  V1  V1 


(C  -  U  >  C  P  +1 


c  . 


VV1 


This  quantity  is  (compare  with  (4.38))  nonnegative  for  1  <  Y  <  3.  Hence  the  minimum 
of  u  1j  —  c.  i,  occurs  for  u_,  •  c m  c  .  We  substitute  this  value  into  (4.37),  and 

must  show: 


.1A"  C1  -  ’A  occurs  for  ui  “  ci 
Jq  ^2  Jq  ^  30 


(4.44) 


(Tt1)(uL  -  777  cL>  >  (y-3)(uR  +  777  cR)  , 


but  we  have  already  shown: 

_2_  XtI  2  Y-1 

,  L  _  ,  R  Y+1 .  R.Y+1  .  L  R.Y+1  ,  R.Y+1 

(4.45)  c  <  (c  )  (u  )  — *>  u  >  (c  )  (u  ) 

and  thus,  by  Holder's  inequality  and  the  fact  that  Y  <  3 


(4.46)  (Y+1)(uL  -  777  cL)  >  (Y+1)(^77  )<cR>1f+1  (“V*1  >  (Y-3)(uR  +  777  cR 


Next  we  note  that: 


(4.47) 


V1  +  Ci0  ,  Y-1  .  . 

VV2 - 5 —  +  V<UV1- V 


cR  Y-1  R  j0  Y-1 

2  +  4  U  +  2  "  4  Ujn 


by  (4.29 )( c) . 


Thus,  using  (4.31)  -  (4.34),  and  (4.38),  we  have: 


(cj  +1  +  cj  1 

■’o  ^0  Y-1  , 

- 2 - —  (Uj0+1  *  V 


3-y  R  2  R  Y+1  2  , 

~  (u  +  — r  c  )  +  —r-  (u.  -  — r  c.  )  , 
4  Y-1  4  30  Y-1  3q 


by  (4.29)(e).  Hence  by  (4.48) 
(4.50) 


<*>  'U3n+V2‘  Cj+V2  "  "  Y-1  do,  Cj  +V  ° 


Thus  the  minimum  of  u.  -  c  .  occurs  at  (c.  ,u  )  =  (c,c),  where 

30+  '2  V  '2  3n 

(4.51)  min( u 


c,  +  1,1  =  0  . 

0 

We  have  thus  verified  that  each  member  of  the  relevant  family  of  sequences  defined  in 


V  v  W 


the  statement  of  Theorem  (4.2)  is  indeed  a  discrete  one  shock. 


!  1 
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An  analogous  verification  works  in  the  discrete  two  shock  case 


We  finally  turn  to  the  proof  of  uniqueness  -  that  these  are  the  only  steady  discrete 
shock  solutions  of  (3.23),  (3.25),  or  (3.26). 

We  begin  with 
Lemma  (4.1) 

Suppose 

/  <3f(w))+dw  +  /  (3f(w))"dw  -  0 

rj  r0+1 

for  j  "  j0  ,  j0  +  1  .  Moreover,  suppose  none  of  the  eigenvalues  of  3f(w^)  either  vanish 
or  change  sign  for  j  -  j„-1,  jQ  ,  jQ+1 ,  jQ+2.  Then  w  -  w  . 

Proof 

Suppose  first  the  eigenvalues  of  3f(w)  are  both  positive.  Then  f(w  )  -  f(w.  ) 

30  V1 

and  by  evaluating  (4.52)  for  j  -  Jg  ,  we  have 

(4.53)  f(w  )  -  f(w  _  )  »  /  3f(w)dw  +  /.  3f(w)dw 

Jq  Jq 

F2  r2 

-  +  -j„)r2<w3  1  +  h  %  "  s  >Vwj  >  -  0  ' 

J0  J0  J0  J0  Jo  J0  j  Jo 

for  some  intermediate  values  w.  ,  w.  ,  and  with  Tr  “|w.  .,-w  ,lr.  , 

"j0  "j0  _j0  V  *  V 1 

Y.  •>  |w  -  w  ijc-  ,  with  c.  ,  c.  bounded  above  and  below  by  positive  constants. 

-'o  ■*  ■'o-  “]0 

This  follows  from  the  mean  value  theorem  for  integrals.  The  vectors  r<w  ),  r  (w.  )  are 

a  b  “  2  _:I0 

of  the  form  (  ^  ),  (  )  for  a,b  >  0  ,  hence  they  are  linearly  independent,  and  it 


follows  that  w 


,  ■  w  ,  i,  ”  w.  •  Thus  if  all  the  eigenvalues  are  positive,  then 

V1  V  '2  In 


If  all  the  eigenvalues  are  negative,  then  an  analogous  argument  gives 

jo  V1  V2 

Next  suppose  the  eigenvalues  are  of  different  signs,  i.e.  -c^  <  , 

1  "  jg  +  r-1,  r  ■  0, 1,2,3.  Then  the  mean  value  theorem  for  integrals  allows  us  to  rewrite 
(4.52)  as 
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(4.54) 


0 


V-i  +  -j)r2<Sj>  +  Ij+l‘Sj+1  ■  -j+1  >ri  (-j+,  ’  ’  0 

for  j  =  jQ  ,  jQ  +  1.  The  linear  independence  of  r^  and  r2  again  Celle  us  that 
lj  -  2j+1  “  0  which  means  wj.  i/2  ’  wj_1'  wj+i  "  wj+  '4  tor  3  “  30  '  j0  +  1  *  ThuB 

30  V  V2  j0+1 

As  a  consequence  of  this  we  have: 

Corollary  (4.1  ) 

Any  steady  discrete  shock  solution  is  eventually  constant  -  i.e.  there  are  indices 
L  R 

J.j,J2  such  that  Wj  =  w  if  j  <  Wj  =  w  if  j  >  J^. 

Now  we  consider  the  one  shock  case.  Let  jQ  be  the  largest  index  so  that  c^  <  u^ 

L 

for  all  j  <  jg  .  This  means  that  w^  =  w  for  j  <  jQ-1  ,  c^  <  u^  and 

V1  ”  Uj0+1  ‘ 

This  also  means  that  /  (3f(w))  dw  »  0  .  Now  if  u  i,  <  c  i.  ,  this  integral 

Jr,  Jn~  *2  'i 


a  nontrivial  linear  combination  of  two  linearly  independent  two  vectors  (by  a  now  familiar 
application  of  the  mean  value  theorem).  Thus  u,  i.  >c,  y,  . 

V  '2  V  * 


This  implies  that  equation  (4.52)  for  j  “  j0  becomes 


(4.55) 


f(w.  )  -  f(wL)  +  /.  ( 3f  (w)  )~dw  -  0  . 


We  wish  to  show  that  u.  - ,  >  cJ  i,  .  Suppose  that  thiB  is  false.  We  use  our 

V  ^  j0  /2 

special  hypothesis  for  the  first  time  -  u.  >  -c  ,  u  .  >  c  .  .  Then 

Jo  1  30  30  '  ■10  ' 

equation  (4.52),  for  j  “  jg  +  1  ,  becomes 

(4.56)  /  +1Of(w))+dw  +  /.  +1  ( 3f  (w)  )”dw  -  0  . 


Our  usual  "linear  independence"  argument,  together  with  the  assumption 
u  i  >  -c  i,  ,  tells  us  that  w^  i ,  “  wu  .  This  is  a  contradiction  since  we  assumed 

j0+Vi  jQ+  V2  jQ+ V2  jo 


u .  ^  y,  <  c,  .y.  and  u  >  c  . 

V  V2  V  V2  30  30 
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5.  ALGORITHMS  AND  NUMERICAL  EXAMPLES 

In  this  section  we  shall  present  results  from  numerical  tests  with  the  one  dimensional 
potential  flow  scheme  developed  in  Section  3.  We  shall,  however,  first  briefly  discuss 
some  other  related  algorithms  in  one  and  two  space  dimensions.  A  forthcoming  computational 
paper  will  contain  experiments  with  some  of  these  algorithms. 

The  one  dimensional  algorithm  which  is  presented  in  this  paper  can  be  used  as  a  basic 
component  in  an  ADI-scheme  for  the  two  dimensional  full  potential  equation.  The 
appropriate  eigenvectors  can  be  derived  analytically  and  the  corresponding  upwind  scheme 
for  the  three  components  can  be  constructed.  Experiments  with  this  type  of  upwind  schemes 
for  the  Euler  equations  are  given  in  [17]. 

All  methods  in  this  paper  and  in  [17]  are  explicit.  Implicit  methods  often  have  the 
advantage  of  being  stable  for  longer  time  steps  (At).  The  linearized  problem  may  be 
unconditionally  stable.  This  is  a  definite  advantage  for  calculations  to  steady  state  and 
also  for  some  transient  problems. 

An  important  part  of  the  computational  work  when  using  an  implicit  method  lies  in  the 
solution  of  linear  system  of  equations.  We  see  two  ways  to  reduce  the  bandwidth  of  these 
systems  and  hence  also  reduce  the  work. 

The  discrete  scheme  for  the  potential  flow  equations  can  be  rewritten  in  scalar  form, 
compare  [8].  There  is  however  no  simple  closed  form  for  the  corresponding  scalar 
difference  scheme  and  approximations  are  needed.  Another  possibility  is  splitting  the 
differencing  and  the  solution  of  linear  systems  into  two  parts;  one  corresponding  to 
forward  and  one  to  backward  differencing.  The  linear  systems  to  be  solved  are  then  block 
bidiagonal  and  the  computational  complexity  is  reduced  compared  to  a  straight  forward  use 
of  Crank-Nicolson  or  fully  implicity  schemes. 

Our  numerical  example  will  be  the  potential  flow  equations  in  one  space  variable  (1.6) 


(5.1) 


[  p  )  ♦  f 

1  u  ;t  v 


pu 


1a  II2  + 


y2  u 


hL. 

Y-1 


Y-1 


0  <  x  < 


t  >  0 


with  initial  and  boundary  conditions 
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(5.2) 


p(x,0)  -  p(x)  ,  u(  x,0 )  -  u(x) 


(5.3) 


p(0,t)  -  PL(t)  ,  u(0,t)  -  uL(t) 

u(1,t)  -  uR(t> 


The  boundary  conditions  (5.3)  in  our  example  correspond  to  supersonic  inflow  at 
x  -  0  and  subsonic  outflow  at  x  ■  1.  Other  combinations  have  also  been  studied  and  they 
give  the  same  qualitative  behavior  except  for  the  fully  supersonic  case.  When  |u|  >  c 
the  scheme  derived  in  Section  3  will  reduce  to  the  standard  upwind  scheme  without 
switches.  Shocks  will  then  not  have  sharp  profiles  as  was  also  the  case  in  scalar 
approximations  [5]. 

The  following  values  we  have  chosen  agree  with  those  in  the  examples  in  [8] . 

Y  -  1.4 

A  -  1/(1. 4  •  1.22)  «  t  i 
PL(t)  5  u^tt)  =  1  . 

The  value  of  u  on  the  right  will  be  changed  between  the  two  examples  presented  here. 
Example  Is  This  is  exactly  the  problem  considered  in  (8]  with  uR( t)  »  0.8.  The  one  shock 
solution  we  will  approximate,  is  characterized  by  the  right  state  PR(t>  -  1.265  and  the 
shock  velocity  VR  -  0.0447.  In  the  scheme  (2.9)  the  above  value  of  pR(t)  was  9lven  as 
numerical  boundary  condition.  See  Section  3  for  the  influence  of  numerical  boundary 
conditions.  See  Figure  5.1  and  5.2  for  computational  results  and  parameters.  The  shock  is 
essentially  resolved  over  two  mesh  points  even  if  the  theory  of  Section  4  does  not  cover 
this  case. 

Example  2;  In  this  test  we  modify  the  outflow  values  in  order  to  achieve  a  steady  shock 

u. ( t )  -  0.728,  p  <t)  -  1.373.  Figures  5.3  and  5.4  display  the  results  of  the 
K  R 

computation.  The  shock  is  resolved  exactly  over  two  mesh  points  and  the  convergence 
history  does  not  produce  any  overshoots. 


2- 


(5.2) 


u(x,0)  »  u(x) 


A*—- 


*■  *  * 


(5.3) 


P(x,0)  -  p(x)  , 

P(0,t)  -  p  (t)  ,  u(0  ,t )  -  U.  (t) 

U  l, 

u(1,t)  -  uR(t) 


The  boundary  conditions  (5.3)  in  our  example  corraspond  to  suparsonic  inflow  at 
x  “  0  and  subsonic  outflow  at  x  «  1.  Other  combinations  have  also  been  studied  and  they 
give  the  same  qualitative  behavior  except  for  the  fully  supersonic  case.  When  |u|  >  c 
the  schesM  derived  in  Section  3  will  reduce  to  the  standard  upwind  scheme  without 
switches.  Shocks  will  then  not  have  sharp  profiles  as  was  also  the  case  in  scalar 
approximations  (5). 

The  following  values  we  have  chosen  agree  with  those  in  the  examples  in  (3) . 

y  ~  1.4 

A  -  1/(1. 4  •  1.22)  -  0.5 
PL(t)  S  S  1  " 

The  value  oi  u  on  the  right  will  be  changed  between  the  two  examples  presented  here. 
Example  1»  This  is  exactly  tne  problem  considered  in  fS]  with  uR(t)  »  0.8.  The  one  shock 
solution  we  will  approximate,  is  characterized  by  the  right  state  PR(t)  »  1.265  and  the 
shock  velocity  Vg  ■  0.0447.  In  the  scheme  (2.9)  the  above  value  of  pR(t)  was  given  as 
numerical  boundary  condition.  See  Section  3  for  the  influence  of  numerical  boundary 
conditions.  See  Figure  5.1  and  5.2  for  computational  results  and  parameters.  The  shock  is 
essentially  resolved  over  two  mesh  points  even  if  the  theory  of  Section  4  does  hot  cover 
this  case. 

Example  2»  In  this  test  we  modify  the  outflow  values  in  order  to  achieve  a  steady  shock 
Ug(t)  ”  0.728,  PR(t)  “  1.373.  Figures  5.3  and  5.4  display  the  results  of  the 
computation.  The  shock  is  resolved  exactly  over  two  mesh  points  and  the  convergence 
history  does  not  produce  any  overshoots. 
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Figure  5.1,  5.2:  Shock  profiles  for  example  1,  section  5. 

Results  after  50  iterations.  Exact  initial 
values.  Courant  number  0.75. 


Shock  profiles  for  example  2,  section  5. 
Results  after  150  iterations,  piecewise 
linear  initial  values.  Courant  number  o.75. 
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